package org.dmikis.task11.engine;

import static java.lang.Math.abs;

import org.dmikis.task11.PhasePoint;

public class RungeKuttEngine extends EulerEngine {
    
    private RungeKuttEngine() {
	super();
    }
    
    public static RungeKuttEngine getInstance() {
	return _engine;
    }
    
    public void run() {
	super._isRunning = true;
	double current_time = 0.0;
	double tmpV1, tmpV2, tmpTetta1, tmpTetta2;
	while (current_time < 50.0) {
	    this._timeStep = 2.0*BasicEngine.TIME_STEP;
	    PhasePoint current_point = points.get(points.size() - 1);
	    do {
		this._timeStep = 0.5 * this._timeStep;
		double m1, m2, m3, m4, k1, k2, k3, k4;
		m1 = this._timeStep*evaluateDtetta(current_point.getV(), current_point.getTetta());
		k1 = this._timeStep*evaluateDv(current_point.getV(), current_point.getTetta());
		m2 = this._timeStep*evaluateDtetta(current_point.getV() + 0.5*m1, current_point.getTetta() + 0.5*k1);
		k2 = this._timeStep*evaluateDv(current_point.getV() + 0.5*m1, current_point.getTetta() + 0.5*k1);
		m3 = this._timeStep*evaluateDtetta(current_point.getV() + 0.5*m2, current_point.getTetta() + 0.5*k2);
		k3 = this._timeStep*evaluateDv(current_point.getV() + 0.5*m2, current_point.getTetta() + 0.5*k2);
		m4 = this._timeStep*evaluateDtetta(current_point.getV() + m3, current_point.getTetta() + k3);
		k4 = this._timeStep*evaluateDv(current_point.getV() + m3,current_point.getTetta() + k3);
		tmpV1 = current_point.getV() + (k1+k2+k3+k4)/6.0;
		tmpTetta1 = current_point.getTetta() + (m1+m2+m3+m4)/6.0;
		m1 = 0.5*this._timeStep*evaluateDtetta(current_point.getV(), current_point.getTetta());
		k1 = 0.5*this._timeStep*evaluateDv(current_point.getV(), current_point.getTetta());
		m2 = 0.5*this._timeStep*evaluateDtetta(current_point.getV() + 0.5*m1, current_point.getTetta() + 0.5*k1);
		k2 = 0.5*this._timeStep*evaluateDv(current_point.getV() + 0.5*m1, current_point.getTetta() + 0.5*k1);
		m3 = 0.5*this._timeStep*evaluateDtetta(current_point.getV() + 0.5*m2, current_point.getTetta() + 0.5*k2);
		k3 = 0.5*this._timeStep*evaluateDv(current_point.getV() + 0.5*m2, current_point.getTetta() + 0.5*k2);
		m4 = 0.5*this._timeStep*evaluateDtetta(current_point.getV() + m3, current_point.getTetta() + k3);
		k4 = 0.5*this._timeStep*evaluateDv(current_point.getV() + m3,current_point.getTetta() + k3);
		tmpV2 = current_point.getV() + (k1+k2+k3+k4)/6.0;
		tmpTetta2 = current_point.getTetta() + (m1+m2+m3+m4)/6.0;
		m1 = 0.5*this._timeStep*evaluateDtetta(current_point.getV(), current_point.getTetta());
		k1 = 0.5*this._timeStep*evaluateDv(current_point.getV(), current_point.getTetta());
		m2 = 0.5*this._timeStep*evaluateDtetta(current_point.getV() + 0.5*m1, current_point.getTetta() + 0.5*k1);
		k2 = 0.5*this._timeStep*evaluateDv(current_point.getV() + 0.5*m1, current_point.getTetta() + 0.5*k1);
		m3 = 0.5*this._timeStep*evaluateDtetta(current_point.getV() + 0.5*m2, current_point.getTetta() + 0.5*k2);
		k3 = 0.5*this._timeStep*evaluateDv(current_point.getV() + 0.5*m2, current_point.getTetta() + 0.5*k2);
		m4 = 0.5*this._timeStep*evaluateDtetta(current_point.getV() + m3, current_point.getTetta() + k3);
		k4 = 0.5*this._timeStep*evaluateDv(current_point.getV() + m3,current_point.getTetta() + k3);
		tmpV2 += (k1+k2+k3+k4)/6.0;
		tmpTetta2 += (m1+m2+m3+m4)/6.0;
	    } while ((abs(tmpV1 - tmpV2) > _error) || (abs(tmpTetta1 - tmpTetta2) > _error));
	    points.add(new PhasePoint(tmpV1, tmpTetta1, current_time += this._timeStep));
	}
	super._isRunning = false;
    }
    
    private static RungeKuttEngine _engine = new RungeKuttEngine();

}
